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The theoretical framework for higher-order correlation functions involving multiple times and mul- 
tiple points in a classical, many-body system developed by Van Zon and Schofield [Phys. Rev. E 65, 
011106 (2002)] is extended here to include tagged particle densities. Such densities have found an 
intriguing application as proposed measures of dynamical heterogeneities in structural glasses. The 
theoretical formalism is based upon projection operator techniques which are used to isolate the 
slow time evolution of dynamical variables by expanding the slowly-evolving component of arbi- 
trary variables in an infinite basis composed of the products of slow variables of the system. The 
resulting formally exact mode-coupling expressions for multiple-point and multiple-time correlation 
functions are made tractable by applying the so-called iV-ordering method. This theory is used 
to derive for moderate densities the leading mode coupling expressions for indicators of relaxation 
type and domain relaxation, which use dynamical niters that lead to multiple-time correlations of 
a tagged particle density. The mode coupling expressions for higher order correlation functions are 
also succesfully tested against simulations of a hard sphere fluid at relatively low density. 
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I. INTRODUCTION 

A complete understanding of the physical processes 
underlying the transition between the high-temperature 
exponential relaxation of density fluctuations of a fluid 
and the non-exponential relaxational profiles observed 
at lower temperatures, especially near the glass tran- 
sition, remains elusive. A number of interesting fea- 
tures have been noted in dense, supercooled systems from 
computer simulation studiesQja, IE El IE IE as weu 
as multi-dimensional NMR|3, Q and video microscopy 
experiments 0, |n| that appear to be related to this 
cross-over from simple exponential to multi or stretched 
exponential relaxation: One typical feature of such sys- 
tems is the appearance of heterogeneously-distributed re- 
gions of the fluid which differ dramatically in their mobil- 
ity and local density. While fluid motions are relatively 
unrestricted in regions of low density, structural rear- 
rangements in regions of high local density at a given time 
have been observed to occur through relatively rapid, 
collective string-like motions0. Furthermore, regions 
which at one time were of relatively low density and in 
which particle motions were primarily fluid-like can be- 
come locally-dense and immobile. 

It is well-known that a variety of different mechanisms 
are consistent with non-exponential relaxation and that 
this relaxation is somehow related to the heterogeneous 
nature of the dense fluid. In one possible scenario of non- 
exponential relaxation, fluctuations of local density relax 
in the same intrinsically non-exponential way, where the 
cooperative motion of particles depends strongly on the 
local environment and is correlated over a long time pe- 
riod or "history" . Another possible mechanism is that 
each region of the fluid rearranges in a less strongly- 



correlated fashion but with different rates, so that the 
observed non-exponentiality is a consequence of the su- 
perposition of different exponential relaxation processes. 
Of course neither scenario may apply over all time scales 
and the mechanism may shift from one that is primarily 
homogeneous to one that is primarily heterogeneous. It 
is equally possible that one type of relaxation may not 
even predominate over another. 

In order to distinguish between these scenarios it is 
helpful to construct quantitative measures which unam- 
biguously signal the presence of a specific mechanism. 
Such constructions can be based on filtersH [H EE El 
that select out sub-ensembles of particles to have specific 
dynamical properties over a sampling period. A simple 
example of such a filter is one that selects out individ- 
ual particles that move either more or less than a critical 
distance over a fixed period of time. One can then ex- 
amine time correlations within these sub-ensembles to 
gain new insight into detailed features of the dynam- 
ics. When time-filters are utilized in this fashion, time 
correlation functions of particles contained in the sub- 
ensemble necessarily involve multiple-time intervals. Fil- 
ters based on single-particle properties then can be ex- 
pressed as multiple-time correlation functions of tagged 
particle densities. 

Given this interesting application of multiple-time cor- 
relation functions, the need for a theory that enables one 
to calculate such quantities from first principles is clear. 
In Ref. |l5| such a theory was developed based on mode- 
coupling theory for multiple-time correlation functions of 
collective densities (i.e., particle number, momentum and 
energy densities). The theory was tested successfully on 
a hard sphere gas at moderate densities in the hydrody- 
namic limit where the relevance of mode-coupling theory 
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is well-established. In this article, the theory is extended 
to include multiple-time correlation functions of arbitrary 
type, encompassing correlation functions involving only 
tagged densities, only collective densities as well as any 
combination of the two types. 

The paper is organized as follows: The mode-coupling 
formalism of multiple-time correlation functions is intro- 
duced in Sec. ^ based on projection operator methods, 
and equations describing the evolution of multiple-point 
densities of arbitrary type are obtained. These equations 
are manipulated to yield expressions for correlation func- 
tions of multiple-time intervals that are local in all time 
arguments. In Sec. IIIII a systematic method of determin- 
ing what types of contributions are the most significant 
for a particular multiple-time interval correlation func- 
tion of tagged particle densities is introduced and the 
leading-order expressions for two and three-time interval 
correlations of a tagged particle density are presented. 
In Sec. IIVI specific applications of multiple-time correla- 
tion functions of tagged particle densities are examined. 
In particular, direct measures of relaxation type and of 
the rate at which solid-like domains become fluidized are 
analyzed and leading order mode-coupling expressions 
for these measures are obtained. Sec. IIVI concludes with 
a numerical comparison of the leading order theoretical 
predictions with direct simulation results for a low den- 
sity, hard sphere system in the hydrodynamic limit and 
Sec. contains a summary. 



II. MODE-COUPLING FORMALISM 
A. Slow variables 

The basic assumption of mode-coupling theory is that 
the long time behavior of all time correlations functions 
can be completely expressed in terms of the evolution of 
a set of slow modes of the system. Although the theory 
doesn't specify the identity of the slow modes, physical 
arguments can often serve as a guide to define a finite set 
of variables that serve as a dynamical basis for the long 
time evolution. For example, since the particle number, 
the momentum and the total energy of a fluid of struc- 
tureless particles are constants of motion, a minimal basis 
set for the slow evolution of the system must include the 
long wavelength modes of densities of these quantities. 

Once the slow basis set has been determined, the slow 
component of an arbitrary dynamical quantity can be 
extracted by finding the projection of the variable onto 
the subspace spanned by the set of slow variables. Simi- 
larly, projecting the variable onto the complement of this 
subspace should yield a fast quantity. 

Projection operators are linear operators, and hence 
only the linear dependence of a dynamical variable on 
the slow variables can be projected out by such a pro- 
cedure. In general, however, one expects the time de- 
pendence of most dynamical variables to be a non-linear 
function of the slow variables. This non-linearity can 



be incorporated into the theory by constructing the slow 
subspace to include the space spanned by all powers of 
the slow variables. In this way an analytic dependence 
of the slowly-evolving component of a dynamical variable 
on the slow variables can be described. The basis of the 
slow subspace is referred to as the multi-linear basis. 

Consider an equilibrium system composed of N par- 
ticles that is described by a translationally invariant 
Hamiltonian TL. The time evolution of any quantity 
(phase space function) C is given by 



C{t) = CC(t), 



(1) 



with £ the Liouvillian operator, which is the Poisson 
bracket with the Hamiltonian. 

To describe the slow evolution of dynamical variables 
within the projection operator formalism, the slow vari- 
ables of the system are taken together in a single vector 
B. As mentioned above, for a structureless fluid, these 
are typically the number density, the momentum density 
and the energy density: 
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Here, r n (t) is the position of particle n at time t, p n (i) 
is its momentum and e n (t) its energy (kinetic and poten- 
tial). It will be convenient to work in Fourier space: 



B(k,t) 



N 



(2) 



where b n [t) = (1, p„(i), e n (t)). 

Since the goal of the mode-coupling theory outlined 
here is to describe the time correlation functions of single 
particles, slow tagged particle densities must be included, 
i.e., 

JVi(r.t) = a(r-n(t)), 

where particle 1 is the tagged particle. In the Fourier rep- 
resentation, the tagged particle number density is simply 



(3) 



The tagged particle density is taken together with B(k) 
into a larger vector A(k). Henceforth, components of A 
will be denoted with a superscript s when referring to sin- 
gle particle densities (which have no summation over par- 
ticles in their definition) and with a superscript c when 
involving collective densities (which have a summation 
over all particles in their definition). 

We require that the slow subspace be spanned by 
all powers of the slow variables. Since the product 
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A r i(k)A r i(q) = Ni(k + q) is just a linear variable, be- 
yond the linear level only products of B are relevant. 
Furthermore, to accommodate correlation functions of 
certain fast variables one might be interested in (e.g. 
Pi(k) = pi0)e 4k - ri (*)), A'(k) is defined as the vector 
composed of all linear slow variables plus any fast vari- 
ables of interest. 

It will prove to be convenient to construct a basis set 
that is orthogonal in the number of factors of the linear 
set A and B, which is guaranteed by the definition: 

Qo = i 

Qi(k) = A'(k) - <A'(k)> = i'(k) 

Q 2 (ki,k 2 ) = i(ki)B(ka) 
1 

- J2(Mki)B(k 2 )Q* Q )*K-}*Q & 

|a|=0 



Here the following notation has been used: 

• (■ ■ • ) denotes the (grand canonical) equilibrium en- 
semble average, which is used to define the inner 
product. 

• A superscript "*" defines complex conjugation. 

• A Greek lower case letter denotes a set of pairs, 
each pair containing a component index and a wave 
vector. 

• \a\ denotes the number of pairs in the set a, the 
so-called mode order. 

• A hatted Greek letter has the same mode order as 
its unhatted variant but is otherwise unrelated. 

• Q a is the same as Qi a i(ki, . . . ,k| Q |). Also, in the 
rest of the paper, the short hand notation of a num- 
ber n instead of a Greek letter will often be used 
to indicate a set of mode order is n. 

• The "*" product involves a summation over the 
pairs of wave vectors and indices in a, divided by 
all ways in which these pairs can be permuted, so 
as to avoid over-counting. Also, for this summa- 
tion to always be well-defined, the wave vectors 
will be summed up to a cut-off k c , thus including 
M = O(N) wave vectors. 

• By definition, K a p = (Q a Q*p), while K^ 1 is the 
inverse of this object with respect to the V prod- 
uct. 

As the definition of Q n contains K mm for all m < n, 
which arc defined using the Q m <n, the above definition 
of the multi-linear basis is really recursive. Note that as 
Knm = if 7i ^ m, the above set is orthogonal in mode 
order: K nrn = (Q n Q m ) — unless n = m. 



We note also that Nx (k = 0) = (exp[ik-ri]— <5ko)| k _ = 
0. This implies that (Nx(0,t)N*(0)) = and, more im- 
portantly, that the elements Q a which have a Nx compo- 
nent with zero wave vector have to be omitted from the 
set to avoid an over-complete basis. 

B. Single time interval correlations 

Let all time correlation functions of the slow (linear) 
variables be taken together into the matrix 

G u (t) = {Qi(t)Qt)*Krf (4) 

Note that Gxx is wave vector dependent, but that this is 
not explicitly denoted here. 

In order to make contact with other versions of mode- 
coupling theory, we first examine the consequences of 
taking only the linear basis Qx into account. The time 
evolution of Q\ is given by Eq. ^ In differential form, 
this means 

Q l {t)=LQ 1 {t) = e a CQ l . (5) 

Defining the linear projection operator 

V X C = {C Ql) * i^i 1 * Qx 

and its complement = 1 — V\ , and using a well known 
operator identity, Eq. can be cast in the form of a 
generalized Langevin equation 

Qx(t)=M 1 E 1 *Qx(t)+ f Mxx{t-T)*Qx{r)dT + (px{t), 
Jo 

(6) 

in which a fluctuating force (f>x{t) = e Vl ct Vi£Qx, a 
static vertex Mf{ = ({CQx} QT) * Kxi an< ^ a Tnemory 
kernel Mxx(t) = ~( l fi(t) ( f*) * ^"ii appear. By taking 
the inner product with Qx, Eq. yields 

Gxx(t)=MZ*Gxx(t)+ f Mxi(t-T)*Gxx(T)dr. (7) 
Jo 

It should be noted that because of translational sym- 
metry, the sum of wave vectors in an average has to add 
up to zero. Thus Gxi involves only one wave vector, 
instead of two. Likewise, the above equation is an equa- 
tion involving just one wave vector. In this sense, there 
is no mode-coupling in the above equation, although, of 
course, the different components of Qx are coupled. 

The memory kernel Mn(i) involves the time correla- 
tion of the fluctuating force ipi, for which the formalism 
as presented thus far provides no method of evaluating. 
But provided projects out all the slow behavior of 
tp(t), (ifx(t)<fil), which contains dynamics orthogonal to 
A only, should be a quickly-decaying function that is well- 
approximated at long times t ^> t m by Mxx(t) ~ 2D5(t), 
for some microscopic time t m on the order of the particle- 
particle collision time. Under such circumstances, Eq. 
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is local in time, i.e., G u (i) = (M^ + D) * G n (t). Un- 
fortunately, in most cases, the correlation function of the 
dissipative force <pi(t) is not a qui ckly -decaying function 

but instead has long time tails[la[l21L3 "^ ue *° ^ ne ^ ac * 
that the linear projection operator does not remove 
the non-linear dependence of tp(t) on the set of slow vari- 
ables B. Hence, one is forced to use the full multi-linear 
basis if equations of motion are to be local in time. 

In the following, Q will denote the vector composed 
of all Q a - As Q is still a phase space function, its time 
evolution is governed by 

Q(t)=CQ(t). 

Using the multi-linear projection operator 

VC ee (CQ*) *K~ 1 *Q, 

(where the "*" now also denotes a sum over mode-orders) 
and its complement = 1 — V, an equation analogous 
to Eq. is found: 

G(t) = M E *G{t) + f M D {t-T)*G{r)dT, (8) 
Jo 

where G(t) is the full propagator defined by 

G(t)=(Q(t)Q*)*K-\ (9) 
the fluctuating force is defined by 

tp(t) = e v " ct V^LQ, (10) 

and the vertices M E and M D {t) are given by 

M E = ({£Q}Q*)*K-\ 
M D (t) = —((p{t) ip*) * K~ x . 

Note that the original goal of evaluating Gu(t) now be- 
comes to calculate the 1 — 1 (or linear-linear) sub-block 
of G(t) while Go a (t) = G a o{t) = S a0 is trivial. 

Now that all powers of B are projected out of the dy- 
namics, one expects p to be truly fast, and its correlation 
function to be approximately a delta function. Defining 
the dissipative vertex as 

/>oo 

M D EE / M°(t)dt, 

Jo 

and the full vertex to be 

M ee M E + M D (11) 

we can write 

G(t) = M * G(t), (12) 

as an approximation to Eq. |H1 

A physical note is in order here: it is assumed that 
there is a separation of time scales, such that there are 



fast correlations which decay on the short microscopic 
scale t m , while the interesting long time behavior occurs 
on slow, 'hydrodynamic' scales th, and we assume t m <C 
th- Eq. E| describes just the slow part, so it is valid only 
after a time (much) larger than t rn with corrections of 
order 0(t m /t h ). 

It is useful to perform a Laplace transform: 

/>oo 

GO) = / e~ zt G{t) dt. 
Jo 

With the initial condition G a p(t = 0) = l a p (the unit 
matrix in infinite dimensions) , Eq. 1121 can be solved for- 
mally in Laplace space, 

G(z) = [z-M]- 1 , 

where the inverse is to be taken with respect to the "*" 
product. By splitting up the matrix M a p into its part 
diagonal in wave vectors M d (i.e., the wave vectors in a 
and (3 are pair- wise equal) and an off-diagonal remainder, 
M° = M - M d , one can write 

G(z) = [z-M d -M°]- 1 = G! 6 (*)*[l-Af°G 6 (*)]- 1 1 (13) 

where the bare, mode-order diagonal propagator G b is 
defined as 

G b (z) ee [z-M d ]-\ 

Because the wave vectors in a and [3 have to pair up in 
M d g, the numbers of wave vectors \a\ and |/3| need to be 

equal as well, i.e., M d and G d are also diagonal in mode- 
order. G b (z) at a particular mode-order n is denoted as 
G b n (z)._ 

The inverse in Eq. 1131 can be expanded to yield 

G(z) = G b (z) + G b {z)*M° *G b (z) 

+ G b {z) * M° * G b (z) * M° * G b (z) + . . (14) 

The 1 — 1 element of G(z) can be written in terms of a 
self-energy that is defined as 

oo 

E(*) EE Y; M °n*G b n {z)*M° nl 
n=2 
oo oo 

+ E E M ?n * G b n {z) * M„» m * G b m (z) * Ml\5) 

7-1=2 771=2 
+ ... 

Note that the summations start at mode-order 2. After 
re-summmg Eq.HSl Gn{z) is related to E(z) by 

G n (z) - [z — Mix — E(z)] -1 . (16) 

This result can be compared to the solution in Laplace 
space of Eq.0 G n {t) = [z - M E - M(z)}- 1 . Thus, the 
self energy is related to the memory kernel by 

M{z) = M D [+T,{z). (17) 
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Since is short lived, the long time tails of the memory 
function are due to mode-coupling effects in E(z). 

It was shown by Schofield and Oppenheim^ij that in 
the thermodynamic limit, the series for the self energy for 
£ can be re-summed, with the result that all bare prop- 
agators in Eq. ^j] are replaced by propagators that are 
completely diagonal in wave vector, but with the restric- 
tion of the sum over intermediate wave vector sets that 
none of them should be equal. In the same paper it was 
shown that that these full diagonal multi-linear propa- 
gators factor in the thermodynamic limit into products 
of full linear-linear propagators. As a result, Eg. 1151 and 
Eg. 1161 combine to a self-consistent equation for Gu(z). 
It should be noted that these results rely on the N- 
ordering technique, which will be discussed in detail in 
section UTTl 



C. Multiple-point correlations 

In Ref. Il5l the functions G a p(t) with either \a\ or |/3| 
bigger than one were called multiple-point correlation 
functions because they can be seen as Fourier transforms 
of correlation function of densities involving more than 
one relative position. Note that if \a\ = \(3\ and the wave 
vectors in a and (3 are fully (pairwise) matched, G a /3 is 
the full propagator at mode-order \a\, and this propaga- 
tor can be written as a product of full-linear propagators 
to an excellent approximation. However, even when the 
wave vectors are not fully matched, G a p is still an inter- 
esting but nontrivial quantity. 

The expression for the multiple-point correlation func- 
tions follows from the general form in Eq. The form 
of the equation is identical to the case of just collective 
densities that was treated in Ref. (section II D). By 
performing the same re-summations, which rely on the 
iV-ordering method to be discussed shortly, one obtains 
the result of Eq. (26) of that paper: 

G a (3 — G aa 'S a '/3 + G aa i * M°,p, * G/3'/3 

+ G aa ' * Ma'S * Gs&' * Mgi pi * G 'pi p 

+ .... (18) 

Here, primed Greek indices have the same wave vectors 
as their unprimed variants, but not necessarily the same 
component indices, i.e., they are fully diagonal in wave 
vector. Thus, Gaa'^a'fi denotes the full diagonal in wave 
vector of the propagator at mode-order \a\. Furthermore, 
there is a restriction on the summation that none of the 
intermediate wave vector sets are equal, which in contrast 
to the notation in Ref. is not denoted explicitly here. 

Eq. I|18|) expresses the multiple-point correlation func- 
tion in terms of the full propagators, which can be ex- 
pressed in terms of the full linear- linear propagators. By 
using iV-ordering (see section IIII|) . one can show that 
contributions in Eo. 1 181 involving G n (z) with n < \a\ and 
n < |/3| are negligible in the thermodynamic limit. 



D. Multiple-time correlations 

The results above consider only correlation functions 
that contain a single time interval. However, the 
case of correlation functions of several time intervals, 
which in general will be denoted (following Ref. IHh by 

r>( n ) 



,ti) or 



Gg } ({ti}) = m Q an (h + - ■ ■+t n ) . . . Qai(tl))*K^ ao , 

is of considerable interest. 

Even very straightforward approaches to multiple- 
time correlation functions lead to expressions involving 
multiple-time correlations of the fluctuating force. Such 
correlation functions are gene rally non-trivial as they are 
not always fast decaying [20j. The essential ingredient 
in deriving mode-coupling expressions for multiple-time 
correlation functions, as remarked in Ref. 0, is that 
" in a correlation function involving fluctuating forces, 
the function decays quickly in a pair of time arguments, 
provided these are well-separated in time." Here, "well- 
separated" means having a time separation much larger 
than the microscopic time scale t m on which (<p(t) <p*) 
decays. For correlation functions of collective densi- 
ties, this argument led to the conclusion that a correla- 
tion function involving several time arguments ti whose 
time evolution arises through the fast evolution opera- 
tors exp[V CU] can be considered fast-decaying in each 
of the times ti provided all times are positive, larger than 
t m , and the evolution operators are applied in succession. 
When these properties hold, all effective times (such as 
t\, t\ + 1%, t\ + t% + t$, etc.) are well-separated and the 
quoted criterion applies. 

On a formal level, the arguments invoked in Ref. [15| to 
derive expressions for multiple-time correlation functions 
of collective densities apply equally well when tagged par- 
ticle densities are included in the slow basis set. Thus, 
by using projection operator techniques, the following re- 
lation can be established: 



Ga n j3{tn — T)Mp an _ 1 s 



o Jo 



(tn-1 ~ r li t n -2, 
fast term in £ n _i, 



where 



Mp aS (r, Ti) =4(Qf 3 QaQl}K7 1 8( T )S(T 1 ) 



drdn 



(19a) 



(foeT c ^V^p{r)Qa)K7l (19b) 



In the limit of long times where n 3> t m , the fast term 
can be neglected and the integrand can be replaced by a 
delta function according to the above criterion, yielding 
an equation that is local in t n , 

G<£ (tn, . . .) = Ga n0 (t n )Mp an _ lS G ( s n ~ 1 \ (*„_!,. . .), 

(20) 
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where M$ a g = dr\ J °° dTMs a e(T,Ti). Different from 
the single time correlations of sections 111 Bl and 111 CI 
where including Qq was of little practical consequence, 
for the multiple-time correlations here it may happen 
that or \S\ is zero. For these special cases, Mf3 a 6 are 
given by Mp a0 = Kp a * and M 0a $ = l aS . 

The recursion relation in Ea. l2()l can be applied as many 
times as necessary to yield a relation between G^ n ' and 
G^\ For instance, for n = 2 and n = 3: 



(3) 

^0/375(^3,^2,^1) = G a C t (t i )MQpgGe rj (t2)M rll \Gxs{ti). 



III. TV-ORDERING SCHEME FOR TAGGED 
AND MIXED CORRELATIONS 

In order to make the infinite series such as in Eq. 1151 
or 111 Eq. E01 tractable, one must develop a systematic 
scheme for analyzing the relative importance of terms 
appearing in the series. In such an approach, the se- 
ries must be analyzed so that simple but accurate ap- 
proximations for the entire series can be formulated in a 
(preferably) controlled fashion. The A-ordering method, 
developed by Oppenheim and co-workers |2ll l22j | as an 
extension of Van Kampen's system size expansion |23f. 
allows such an approach for correlation functions on the 
hydrodynamic length scale in systems of moderate den- 
sity removed from any critical point. In the A-ordering 
approach, one assigns a factor of N (the number or aver- 
age number of particles) to any cumulant of multi-linear 
densities based on the assumption that the each cumulant 
containing n linear densities scales with the system size 
as A(£/ 'a)^™" 1 -*, where £ is the correlation length of the 
system and a is the average distance between particles. 
As shown in Ref. flfi, the A-order of an arbitrary correla- 
tion function of basis functions Kaa — (Q&Qa) depends 
on the nature of the densities and the number of matched 
wave vectors in the sets a and a. The subscripts like 



{Ax(k-qi- 



-Q\a\-i),Bi(qi) 



H 



denote sets of wave vector-dependent densities, where 
each Bi is a component of the set of collective slow vari- 
able defined in Eq. and Ai is a component of the ex- 
tended linear basis set A. If the component A\ in the set 
a is a single particle density, we denote the corresponding 
multi- linear density by a superscript s, Q s a ; if this compo- 
nent is a collective density, then the dynamical variable is 
represented by Q a - In Ref. flfi it was demonstrated that 
the leading A-order terms arise from matching as many 
wave vectors as possible in the sets a and d, yielding the 
following estimates for A QQ , 



K c - C 

OLOL 

K s - C ~ 

act 



ATM 
A^H-i 



K c - S 

a a 

Kit 



(22) 



where, for example, the superscript "ss" denotes that 
both sets a and a contain a single particle density. More 



importantly, it was also demonstrated that the iV-order 
°f ~ -/Vl Q l -1 only when the single particle densities 
are in the same matched set, and the A-order drops by a 
factor of N when the wave vectors of the single particle 
densities do not match. 

For the inverse of K, the N ordering in Ea. 1221 implies 



AT-M 



^ AM"! 
jy-M+l_ 



(23) 



When considering the N order of an expression that 
(21) contains the "*" product, one needs to separately con- 



sider its order in M (the number of wave vectors summed 
over, see the seventh point at the end of section Hi Ajl . 
That is, when the "*" product is between multi-indices 
of order n, one has in principle n — 1 sums over wave 
vectors and thus an effective factor of M n ~ x (one loses 
one sum over a wave vector because by translation sym- 
metry all wave vectors have to add up to zero inside an 
average). Although M is of order N, these orders of M 
are counted separately from the N ordering, for the fol- 
lowing good reason: For fluids of moderate density only a 
small fraction of the wave vectors in the sums really con- 
tributes significantly. Therefore, rather than taking M 
as counting the precise number of wave vectors summed 
over, it makes more sense to take M as the number of 
wave vectors that contribute substantially to the sum. 
This results in a small value of M [which is nonethe- 
less still O(N)}. The 0(1) parameter M/N is called the 
mode-coupling parameter. As stated, it is typically small 
(e.g. 10~ 5 ) for fluids of moderate density and away from 
any critical point. 

To illustrate this, consider for instance K^*K^, which 
is seen to have an A-ordering of 0(A 4 ). There is a pos- 
sible summation over M wave vectors in * , but 
in fact the leading order estimate of 0(A 4 ) requires a 
matching of wave vectors which gets rid of the factor M. 
As a result, the leading term of * A||, is just 0(A 4 ), 
with (possible) correction terms of order 0(MN 3 ), which 
compared to the leading term are of relative order M/N, 
i.e., of the mode coupling parameter. 

In general, for any expression one can determine 
the values of m and n such that it is of order 
0(M m N n )=0((M/N) m N n - m ). If m > n, such an ex- 
pression vanishes in the thermodynamic limit A — * 00. 
Below, when adding expressions of different M and A 
order, expressions which vanish (relatively) in the ther- 
modynamic limit will be omitted. The remaining terms 
with the smallest power of M/N (typical (M/N)°) will be 
referred to as the leading A -order terms and the terms 
of one power higher in M/N will be called the leading 
correction terms or the first mode-coupling corrections. 



A. TV-ordering of single time interval correlation 
functions 

Using these principles, it was shown 1^ that the A- 
order of the (multi-linear) vertices in Eq. 1111 (in their 



7 



maximally matched form) are: 



M 
M 
AT 



so 

cs 
80 
cc 
50 



0{N-^- & )) 
O(N ) 



if 8 < f3 
if S>(3 



M\l = N- l O{Ml%) 



Following these lines of analysis, the matrix M in Eg. 1111 
and the normalized single-time interval correlation func- 



tions in Eq. 0are[Tf 

J™dt{ip a (t) <p*p) 



in terms of N a/3 — ({£ Q a } Q*, 



M c a c = N^Kf-^l + OiN- 1 )] 



00 



[z-M cc {z)]- a l 



(24) 



from which we see that single particle modes do not con- 
tribute to the dynamics of single time interval correlation 
functions of the collective modes. Similarly, since the su- 
perscripts ss in M ss in Eq. [21 only mean that one of 
the components in the each of the indices of M ss has 
a single particle character while the others are collec- 
tive, the collective modes do influence the dynamics of 
correlation functions of tagged (single) particle densities 
through mode-coupling corrections to transport coeffi- 
cients. For instance, the generalized self-diffusion coeffi- 
cient -D(k, t) is renormalized through the terms in Ea.1151 
with the bare propagator replaced by the full one (see 
the discussion below Eg. I17fl . i.e., 

D(k,t) =D B (k, t)+£ ss (k,t) 
S ss (k, t) = J2 M 5 (JVi (k) ; JVi (k - q), Bi (q) ) F s (k - q; t) 

x Gl\ {Bi (q) ; Bj (q) ; t) M|; (jV x (k - q) , B s (q) ; TV X (k)) 



where D B = is the "bare" diffusion coefficient with 
weak k and t dependence and F s (k, t) = (TV^k, i)TV*(k)) 
is the self-part of the dynamic structure factor. Care- 
ful analysis of the mode-coupling contributions to tagged 
particle transport coefficients is essential in the deriva- 
tion of observed relations between quantities such as the 
self-diffusion coefficient of a Brownian particle and the 
viscosity of the fluid [H El El ■ 



B. TV-ordering of higher-order vertices and 
correlation functions 



The analysis of the TV-ordering of higher-order vertices 
involving mixed tagged and collective particle densities 
follows by induction as outlined in Ref. Using the 
TV-ordering results for the single time interval correla- 
tion functions and the relation between the single-time 



interval correlation functions and the multiple-time cor- 
relation functions in Eq. [21 one obtains the following 
TV-order in the maximally matched form of the higher or- 
der vertices in which the central index is a linear tagged 
density: 



M css 

iCrsss 
M -yl8 



0(TV-( 5 -x>) 

OiN- 1 ) 
O{N ) 

OiN- 1 ) 
O{N ) 

{ O(TV ) 

o(m;i). 



if 7 < 6 
if 7 = 8 
if 7 > 5 

if 7 < 8 
if 7 = 8 
if 7 > 8 

if 7 < 8 
if 7> 6 



Furthermore, the higher order propagators G 7 is obey 
the same TV-ordering rules as the higher-order vertices, 
namely 



G lU = 0(M 7 i 4 ) 



(25) 



Using these results, we see that the two-time interval 
tagged particle correlation function (with time intervals 
t\ and i 2 ) reduces to a particularly simple form to leading 
TV-order: 

G&te.ti) = <g 1 (k-q,< 1 + t 2 )g s 1 (q,i 1 )g s 1 (k)*) 

= G^(t 2 )M 1 s ffGf 1 (t 1 ) + 0(TV- 1 ), (26) 
with leading order correction terms 

GS(i a )M|f2G5'(*i) + GtK^MiltGiKh) 

+Gl s 1 (t 2 )M^G s 2 s 1 (ti) (27) 

where Gf 2 (t) and G 2 f(i) are higher-order, single time 
interval correlation functions. Note that in Eq. [2U and 
Eq. |2I the matrix indices, as more fully indicated in 
Eq. 1211 have been suppressed for notational simplicity. 
Given the definition of Mi 12 in Eq. ED one sees that 
only the second, dissipative term in Ea . I19bl contributes, 
which is O(fcg)- Thus the first term in Eq.|2|is, in orders 
of fco, the leading correction term, while the others are 
0{kl). 

For the three-time interval correlation function, the 
leading TV-order contribution is 

Gnnfe, h, h) — Gll(t3)Mli Goo(t2)M l{GH(ti) 
+ Gil (h)M^(GH (t 2 )MHfGtl (h), (28) 

with the following three terms contributing at the order 
of the first mode-coupling corrections: 

<?f 2 (ts)MitlG s 2 l (t 2 )MlHGll fa ) 
+ G" (t 3 )MSiGll(t a )MffiG% (h) 
+ GS(t a )M 2 'SG35(ia)MS5G5;(ti) + O(k 2 ). 



(29) 



The 0(ky) correction here in fact consists of seven 
more correction terms, all variants of the form 
Gi2(t3)M2iiGu(t 2 )MniGii(ti) and all 0(*g). 
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IV. APPLICATIONS OF HIGHER ORDER 
CORRELATION FUNCTIONS 

It is often difficult using single time-interval correla- 
tion functions to discern what features of the underly- 
ing collective dynamics of slowly relaxing systems leads 
to qualitative features in the relaxation profile in glassy 
and frustrated systems. Typically, single time-interval 
correlation functions in frustrated systems display non- 
exponential time decay, often exhibiting a two-step re- 
laxation processes associated with caging effects and co- 
operative flow through heterogeneous dynamics 0, 0, @ 
on long time scales. Given the heterogeneous nature of 
the dynamics in such systems, it is natural to ask what 
type of relaxation processes lead to this non-exponential 
time signature. 

In a series of articles 0, 0, 0, Heuer and co- 
workers have examined the information content of higher- 
order correlation functions to assess how detailed features 
of the dynamics correspond to aspects of the relaxation 
in glassy systems. In particular, Heuer et al. have fo- 
cused on multiple-time correlation functions designed to 
probe relaxation typefy ^| as well as rate memor»|12] 
associated with the persistence of slow particle motion in 
supercooled liquid systems. The basic idea of the multi- 
ple time correlation approach is to examine correlations 
of particle motion over several time intervals, separat- 
ing out distance and directional correlation [l3j. In this 
fashion, one effectively devises time filters that extract a 
particular feature of the dynamics to be analyzed. The 
fundamental building block of the time correlation func- 
tions is the (real part of) tagged particle density at time 
interval Atoi = t\ — to defined to be 

/(t ,*i) = cos(k-(r(ti)-r(t ))) = cos(k-Ar 01 ) 

= Re(iVi(k,t 1 )iV 1 *(k,to)) (30) 

whose ensemble average gives the incoherent scattering 
function F s (k, Atoi ) ■ Intuitively, F s (k, At) measures the 
fraction of particles moving a distance less than 2ir/k over 
the time interval At, and hence /(to, ti) can be viewed as 
a time-filter selecting out slowly-moving particles. The 
essential idea in identifying the relaxation type is to con- 
sider how the motion of slow-particles is correlated over 
subsequent time intervals. In the purely heterogeneous 
scenario, one expects that motion in subsequent time in- 
tervals has no direction dependence (no back-and-forth 
motion). On the other hand, for the purely homogeneous 
scenario, one rules out a distance dependence in subse- 
quent time intervals to exclude the presence of different 
mobilities. In order to characterize these limits, it is help- 
ful to define the three-time correlation function fl3| 

F 3 (At 01 ,Ati 2 ) = (/(t ,ti)/(ti,t 2 )) 

= (cos(k- Ar i)cos(k- An 2 )).(31) 

In the homogeneous limit, the projected distance k- Ari 2 
along k is independent of k • Argi, and hence the three- 



time correlation function factors to 

F 3 (Aioi,Ati 2 ) = (/(to, *i)> * 2 )) 

= F s (k,At 01 )F s (k : At 12 ), (32) 

which suggests defining an indicator function for homo- 
geneous dynamics 

if oro (k, At 01 , At ia ) = F 3 (At 01 ,At 12 ) 

-F s (k,At Q1 )F s (k,At 12 ) (33) 

that vanishes in the homogeneous limit. On the other 
hand, in the case of purely heterogeneous dynamics, con- 
sider the indicator|13j 

if et (k, Atoi, Ati 2 ) =F s (k, At w + At 12 ) 

- F 3 (At Q1 , At 12 ) (34) 
= - (sin (k • Ar 01 ) sin (k ■ Ar 01 )). 

Since the direction of the motion in subsequent time in- 
tervals is not correlated in the heterogeneous limit, the 
right hand side of Eg . I3*H vanishes . Note that both indica- 
tor functions make use of the three-time correlation func- 
tion of the tagged particle density and can be expressed 
in terms the multiple-time propagators of the previous 
section as 

F 3 (k,At 01 ,At 12 ) = iF s (k,A% + Ati 2 ) (35) 

+ i (iVi (k, t +h+ t 2 )N x (-2k, t + h) Ni (k, t )*) , 

where the second term on the right hand side is a spe- 
cial case of the more general propagator in Eq. 127)1 
defined with q = 2k and t = and the tagged 
particle densities corresponding to the number density. 
Both measures of relaxation type have been success- 
fully tested 26] on simple I-dimensional model systems in 
which the dynamical rules governing motion of a particle 
are constructed to be inherently heterogeneous (an en- 
semble of particles each moving with constant but differ- 
ent jump rates) or homogeneous (a collection of particles 
in which particles hop between sites with two different 
site-dependent rates). 

As the function /(to, ti) acts as a slow-dynamics filter, 
it can be used as a means of selecting a sub-ensemble of 
the full system. As an alternative to F 3 om and F 3 et in 
studying the heterogeneous nature of the dynamics, one 
can then examine examine how long particles that are 
initially in the slow-dynamics ensemble remain in this 
ensemble to get an idea of how long solid-like domains 
persist in supercooled and glassy systems. Such filters 
are also useful to try to rigorously map deterministic sys- 
tems onto simplified models of glassy behavior, such as 
facilitated spin models^ IS El Ij^ljjlliS- A suitable 
measure of the lifetime of solid-like domains can be de- 
fined by constructing the four-time correlation function 

C^XhMMM) = (f(to,h)f(t 2 ,t 3 )) 

= (cos (k • Ar i ) cos (k • Ar 23 )) . 
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Generally, it is sufficient to look at a time filter over a 
fixed period At = t\ — t = £3 — t 2 , where the waiting in- 
terval <2 = ti + r between subsequent applications of the 
time filter is r. For large waiting times r, once expects 
that only a random selection of the particles initially in 
the slow-ensemble remain in the slow ensemble so that 
(7(4) _ > F s (k, At) 2 as r — > 00. It therefore is logical to 
focus on the fluctuation of f(to,ti) defined as 

C (4) (Ai,r) =(/(0, At)f(r + At, t + 2At)) — F s (k, At) 2 . 

If At is chosen to be shorter than the inverse of the (typ- 
ical) relaxation rate of solid-like domains then the time 
scale r at which this decays to zero can be interpreted as 
the domain relaxation time. 

(3) 

Note that this quantity is related to the G\^ llt ex- 
pressed in Eq. [2H1 and Eq. via 



C^(At, 



1 r 

"4 



_G { l Kli (At,r,At)+G^, KK (At,r,At) 

+ G^, K , R , (At, r, At) + G^ KK , K , (At, r, At) 
-F s (k,At) 2 , (36) 
where k = {Ni(k)} and k' = {JVi(-k)}. 



A. Calculation of domain relaxation rate and 
relaxation type indicators 

Application of the mode-coupling theory of multiple- 
time correlation functions developed in Sec.[n]to evaluate 
the domain relaxation rate via Eq. 1361 or the relaxation 
type indicators defined in Ea. l33l and Ea. l34l reauires com- 
plete specification of the slow basis set variables. As the 
indicators are of significant interest in dense supercooled 
liquids in which F s (k, t) exhibits non-exponential decay 
on molecular length scales, the relevant slow modes must 
describe the long-time evolution of density fluctuations 
for wave vectors k near the peak in the static structure 
factor. Clearly the dynamics at such short length scales 
is outside the regime of hydrodynamic theory for which 
one has a good idea of what constitutes the slow modes 
of the system. For dense systems, however, there is solid 
evidence from the theory of hard sphere liquids[33|, 03 
of the existence of short-wavelength "collective" modes 
that are significantly slower than other "kinetic" modes 
of the system. These collective modes are generaliza- 
tions of the hydrodynamic tagged particle and heat den- 
sity modes to finite wave vectors. The application of the 
mode-coupling theory outlined here to molecular length 
scales is challenging due to the difficulty in evaluating 
the contribution of the fluctuating forces (p a (t) to the 
coupling vertices M Q ^ 7 , and requires new input from ei- 
ther kinetic theory or simulation. Work along these lines 
is in progress. 

In order to get a feeling of what the mode-coupling 
predictions of the correlation functions defined above 
look like, we focus on a moderately dense system (in 



fact relatively dilute compared to a glass) and exam- 
ine these functions in the hydrodynamic limit, as was 
done in Ref. |35[ For such a system, it is sufficient to let 
the set of slow modes be composed of the tagged par- 
ticle number density fluctuations A^i(k) and the collec- 
tive hydrodynamic densities, namely the number density 
fluctuations N(k), the longitudinal momentum density 
Pi(k) — k ■ P(k), and the orthogonalized energy den- 
sity fluctuations H(k) (see Ref. |35| for the precise defi- 
nitions of these variables for a hard sphere system). For 
this specific choice of basis set, the time-derivative of 
the tagged particle number density does have a fluctu- 
ating component ip s N (k, t) that contributes to the M[(( 

vertex. However since the time derivative of N\(k) is 
proportional to k = |k|, one expects these "dissipative" 
contributions to be relatively unimportant in the hydro- 
dynamic limit compared to the non-dissipative couplings 
(QaQpQf) * Kgg~ ' ^ ne higher-order vertices necessary 
to calculate the multiple-time correlation functions for 
the indicator functions to leading A^-order are therefore 



MUl = <A\(k ~ q)A>i(q)A>r(k)) + 0(k 2 ) 
= 1 + (k 2 ) 



(37) 



(Q|(k q qi,qi)A\(q)Qr (k - q'/, q'/)) 

^ll-^k-q^k-qUi) 
+ O{k 2 ) 

^(qi)*^n _1 (qi)<Sqi^ 



O (k 2 ) , 



(38) 



where we have used the factorization properties [l9| of 
multiple-point correlation functions. In the equations 
above, kg represents the largest wave vector present in 
the correlation function. The static part of the multiple- 
point vertices an d M-ffJ coming from the time in- 
tegral of Eq.Eil vanish since Qf(k)Qf(ki) = Qf(k+k x ) 
and (Q2Q1) — by construction. The lowest-order con- 
tribution in wave vector to these vertices therefore comes 
from the time integral of Eq. I19bl and is O(fco)- 



1. Relaxation type 



Combining these results with the leading A^-order ex- 

(2) 

pansion terms of G\H and insertion into the expression 
for F3 in Eq. |23 yields 



F 3 (k,t u t2) =-F s (k,At i+At u ) 



F s (k, At 12 )F s (k, At i) 



+ -F™ c (k,A%,Ai 12 ) 



(39) 



where F™ c is the first mode-coupling contribution to F3 . 
Note that from Eq. 1271 we see the mode-coupling correc- 
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tions involve the evaluation of terms such as 
X)GS (fa (k) ; fa (k - q)B 4 (q) ; At ia ) 

i,q 

* G£ (JVx (k - q)Bi (q) ; fa (k) ; At i) ■ (40) 

In Eq. 1401 the sum extends over the three hydrody- 
namic collective variables -Bj(q) and Gfl, G|f denote the 
multiple-point mixed tagged/collective correlation func- 
tion 

Gr 2 (N 1 (k);N 1 (k- q )B l ( q );t) 

= (jV 1 (k,i)Q5 !jVlfl '*(k-q > q))/<B < (q)B?(q)) 

provided the collective densities are orthogonal 
<Bi(k)S;(k)) = %(Si(k)S|(k)). Similarly, G 2 i is 
defined as 

Gf^iCk-q^q^iViCk);*) 

= (Q^ lBi (k-q,q;t)iVr(k)). (41) 

Using the mode-order expansion Eq. the multiple- 
point correlation function G21 can be approximately 
written as 

G^ Bi;Ari (k-q,q;t)= / dr F s (k - q,t - r) (42) 
Jo 

x £ G Bi Bj (q, t - r)7< lB ' ;Wl (k - q, q; k)F s (k, r), 

with a similar expression for G12. Note that in Ea. 1421 
there is an implicit sum over collective mode index j. In 
practice, it is often convenient to work in a basis set in 
which the matrix of collective linear-linear (normalized) 
correlation functions is diagonal. In the hydrodynamic 
regime, this corresponds to choosing the Bi set to be 
composed of the hydrodynamic sound and heat modes. 

Returning to the problem of calculating relaxation 
type, inserting Eq. 1391 into the definition of the indica- 
tor functions yields 

F*° m (k, At 01 , Ati 2 ) = i[AF 2 (A%, At 12 ) + F™ c ] (43) 
Ft\k, Atou Atia) - ^[AF 2 (At 01 , At 12 ) ~ F™ c ] (44) 
where 

AF 2 (At i,Ati 2 ) =F a (k,Ai i + Atia) (45) 
-F s (k, At i)F s (k, At i2 ). 

In the hydrodynamic limit and ignoring mode-coupling 
effects, one expects that an exponential decay of the in- 
termediate scattering function of the form F s (k, t) ~ 
exp(— D|k| 2 t), where D is the self-diffusion coefficient. In 
this case, AF 2 = and both indicators are approximately 
zero. On the other hand, if F s (k,t) has a non-trivial 
relaxation profile, perhaps of the form of a stretched 



exponential F s (k,t) ~ exp{ — [t/r(k)}^} where (3 is the 
stretching exponent, then AF 2 7^ 0. It is therefore evi- 
dent that mode-coupling effects are absolutely essential 
to distinguish between homogeneous and heterogeneous 
types of non-exponential relaxation. Although this re- 
sult is obvious since mode-coupling must be invoked to 
give rise to non-exponential relaxation in the first place, 
it is the difference in the time dependence of AF 2 and 
the mode-coupling corrections F™ c that allows one to 
distinguish between the two relaxation types. On the 
molecular scale, one also anticipates a contribution to 
these expressions from the dissipative part of Mm that 
corrects the simple factorization result at lowest iV-order 

(iVi(k-q,t 1 +t 2 )iVi(q,ti)iVr(k)) 

»F,(k-q,i 2 )F,(k,t 1 ) [l + 0(*g)] . (46) 

Note that the effect of these corrections is not to modify 
the time behavior but to modify the wave vector depen- 
dence on the right-hand side of Eq.0Tj] This modification 
only makes sense when t\ and ta are large compared to 
the microscopic relaxation time t m corresponding to the 
time scale after which the instantaneous approximation 
implicit in Eq. [2U] is valid. 

2. Domain relaxation 

The rate of domain relaxation can be similarly cal- 
culated using the iV-ordering scheme to simplify the 
multiple-time correlation functions in Eq. 1361 Using the 
form of the higher-order vertices in the hydrodynamic 
limit and G 00 = 1, GH(fa(k = 0); fa(k = 0); r) = 0, 
and Gff (fa(k); fa(k); At) = F s (k,At), one obtains 

G^ KK (At, r, At) a F s (k, At)F s (k, At). (47) 

The same holds for the other four time correlations in 
Eq. therefore the leading orders of the multiple-time 
correlations are canceled by the last term in Eq. [2U so 
that to leading order G^(At, r) is zero. The first non- 
zero contribution arises at first mode-coupling order, as 
given by Eq. 1291 However, the terms in that equation 
that involve Gai(-/Vi(q), Bi(— q); iVi(0); f) are zero be- 
cause iVi(0) = 0. Thus, only the last term from Eq. 1291 
survives. Taking together the contributions from the four 
higher-order terms in Eq. |2U and using that G|| factor- 
izes to leading order, gives the mode coupling result 

(At, r) = Re(j2 G ^ & ( k ) i & ( k " ^ '> At ) 
xf s (q,r)G u '^(q,r) 
x G^ifaik-cOBj^N^At)^ 
+ O(k 2 ), (48) 
where Eq. 021 ma y be used for Gf 2 and G|f . 
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As this result already somewhat suggests, for dense, 
supercooled liquids, one anticipates that the conversion 
time of solid-like domains to liquid-like domains of typi- 
cal length scale I ~ 2ir/k also occurs on time scales corre- 
sponding to the wave vector-dependent a-relaxation time 
of F s (k,t). 



B. Numerical analysis of multiple-time correlation 
functions of tagged densities in a hard sphere system 

In order to validate the mode-coupling theory approach 
to multiple-time correlation functions of collective densi- 
ties, extensive simulations of a hard sphere system at 
moderate densities were carried in Ref. |35l It is straight- 
forward to extend these simulations to incorporate calcu- 
lations of tagged particle densities to test the simple fac- 
torization result for multiple-time correlation functions 
of the tagged particle number density in Eq. I4f>l and to 
examine how well the theory predicts higher-point corre- 
lation functions of mixed tagged/collective densities that 
are necessary to compute the mode-coupling corrections 
to the leading jV-order factorization. All simulation re- 
sults presented in this section were obtained using the 
simulation methodology described in detail in Ref. |35j on 
a hard sphere system of a relatively low reduced density 
p* = p/p c = 0.1, where p c is the density at close pack- 
ing. The size of the periodic system was chosen to have 
cubic box lengths L x = L y = L y = 47.3361 such that 
the simulation system contained N = 15000 hard-sphere 
particles of mass m = 1 and diameter a = 1 at the chosen 
density. For this system size, the smallest dimensionless 
wave vector kga = 2™/ ' L x = 0.132736 so that all quan- 
tities examined are roughly in the hydrodynamic regime. 
For this system, the mean collision time calculated from 
Enskog theory [36fl is approximately t e = 1.42417 at an 
inverse temperature (3 = 3, while the mean- free path 
is l e = 1.85561 so that k l e = 0.246305. Under these 
conditions, the estimated relaxation time of F s (ko,t) is 
r(fe ) ~ (Dkl)- 1 « 55i e , where D = 0.725586 is the self- 
diffusion coefficient calculated from the Enskog theory 
result 



D = 



3w, 



8v J)nmg(a)a 2 



(49) 



where w s = 1.01896 and g{a) is the radial distribution 
function at contact. For a hard-sphere system, g(a) can 
be estimated using the Carnahan-Starling equation of 
state|3_jj and the expression for the pressure p of a hard- 
sphere system 



1 + T] + g 2 — rf 



p (1-77)3 



= 1 + bpg(a), 



(50) 



where b = 2ira 3 /3 and fj = irpa 3 /6. The simulations 
were run for at total time of approximately 1200 r(fc n ) 
to insure reasonably good statistics. Following Ref. 139 . 



statistical uncertainties were estimated using the sym- 
metry properties of the correlation functions. In this ap- 
proach, the statistical uncertainty for a real correlation 
function was constructed from a histogram of the values 
of the imaginary part of the complex correlation function, 
which vanishes on average, to determine the 96% confi- 
dence intervals. To simplify the comparison between the- 
oretical predictions and the simulation results, all wave 
vectors were taken to be co-linear so that k ■ q = kq. 
To further improve statistics, the wave vectors were in- 
dependently taken along the three principal directions x, 
y and z of the cubic simulation box and averaged. 



1. Two time-interval correlation function 

One of the central results of Sec. IIVI the simple fac- 
torization of multiple-time correlation functions of the 
tagged particle number density (see Eq. I46fl is simple to 
verify numerically. To validate this prediction, numeri- 
cal calculations of the correlation function g( 2 >{t\,t2) = 
(Ni(k— q, ii+i2)Ai(q, ti)N^(k)) were computed as func- 
tion of wave vector combinations (k, q) = (to, n)k for to 
and n values ranging from 1 to 3 with m n. To sim- 
plify the comparisons, three different time cross-sections 
(ti, £2) for each pair (to, n) were calculated, namely (t, t), 
(t,3t) and (3t, t). These particular choices of t\ and £2 
are the simplest to implement when simulation data is 
stored in standard linear array data structures. The 
results, shown in Fig. ^ show that the product g^p = 

F s (k— q, r; 2 )-F s (k, t\) approximates </ 2 '(ii,i2) very well 
for all choices of wave vector combinations and all time 
cross-sections. However, given that the statistical uncer- 
tainties for both g( 2 ^ and g^p are on the order of 0.002, 
it is clear that the numerical results indicate a small but 
systematic difference between these two quantities which 
reaches a maximum a short times t ~ 5i e , as can be seen 
in the panels on the right hand side of Fig. ^ 

This difference should be well-reproduced by the mode- 
coupling correction term in Eq. I4UI although this has not 
yet been verified. 



2. Multiple-point correlation function 

Although it is certainly possible to calculate the lead- 
ing order mode-coupling correction term in Eq. I4UI an- 
alytically using hydrodynamic forms for the collective 
correlation functions Gf{(t), the calculation is tedious 
due to the sum over different collective modes -B,;. How- 
ever in order to validate the applicability of the mode- 
coupling theory in the calculation of the first order correc- 
tion terms F"™ c , we examine the quality of the predicted 
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form of the multiple-point correlation function 

Gri(iV 1 (k-q)JV(q);iVi(k);t) (51) 
= (iV x (k - q, t)N(q, t)N* (k) ) - |^F S (k, i), 

where 5(q) = (N(q)N* (q)) is the static structure factor, 
against direct numerical calculation (note: in Ea. 1511 we 
took the thermodynamic limit). From Eq.021 we see that 
this correlation function is approximately given by 

G^ N ^(k- q ,q;t)= f : drF.(k-q,t-T) (52) 
Jo 

x (q, t - r)M^ B ^ (k - q, q; k)F s (k, r), 

where Bj runs over the set {N, P;, 77} or some variant of 
it. Note that with this choice of collective densities, the 
computation of Eq. 1521 requires the calculation of 3 cou- 
pling vertices M% N ' Nl , M^ P ^ , and M^ H ^ , as well 
as the input of the collective linear- linear time correlation 
functions G^ , Gjj and G^ ff . In principle, analyti- 
cal forms of the Gn time correlation functions can be 
utilized in the hydrodynamic regime with transport coef- 
ficients either fitted from simulation data or taken from 
kinetic theory. Given the simplicity and ease of numeri- 
cally calculating the Gn with excellent precision, we have 




10 20 30 40 50 10 20 30 40 50 



Time 

FIG. 1: A comparison of the multiple-time correlation func- 
tion g^ 2 \ the factored form grj 2 ' and their difference as a func- 
tion of wave vector pair (m,n) and time cross-section (ti,ta). 
The left hand side panels contain the simulation results for 
</ 2 ' (unconnected symbols) and g^ (lines), and the right 

hand side panels contain plots of — g^ . In all panels, the 
unconnected dots, crosses and triangles correspond to the sim- 
ulation results for time cross sections (t,t), (3t,t) and (t,3t), 
respectively. The results in the top, middle and bottom rows 
are for wave vector sets (1,2), (1,3) and (2,1), respectively. 
For clarity in the figures, the 96% confidence intervals, esti- 
mated to be roughly 0.002, have been omitted. 



chosen to evaluate the convolution integrals in Ea.l52lbv 
numerically integrating simulation data using a version 
of Simpson's rule that allows interpolation of data. From 
time-inversion symmetry, = for Bj = N or H, and 
it is evident that the only coupling at "Euler" order arises 
for Bj = Pi. However, as noted for the case of multiple- 
point correlation functions of purely collective densities 
in Ref. |35j, the inclusion of the additional couplings to N 
and H arising at "dissipative" order is essential if 
quantitatively accurate predictions are desired. It may 
appear at first glance that the dissipative contributions 
are negligible in the limit of small wave vectors since they 
introduce additional factors of k^. In fact the overall or- 
der of the multiple-point correlation functions is deter- 
mined by a wave vector-dependent prefactor multiplied 
by the convolution of the two-point, single time interval 
correlation functions G\\. The time convolution of these 
correlation functions can give rise to additional factors 
of wave vector depending on their symmetry properties, 
so that the overall contribution of the Pi and N or H 
coupling are comparable in magnitude. 

The evaluation of the dissipative contribution to cou- 
pling vertices AI^ lN ' Nl and M^ lH ' Nl requires external 
input. In the appendix, these vertices are evaluated in 
the small wave vector limit by relating the dissipative ver- 
tices to the Enskog self-diffusion coefficient and its deriva- 
tives with respect to thermodynamic parameters. The 
predictions are therefore free of any adjustable parame- 
ters and constitute a rigorous test of the mode-coupling 
theory. The results of the comparison are shown in Fig. EI 

Note that although the data are a bit noisy, the theo- 
retical predictions generally fall within the confidence in- 
tervals of the simulated data except at very short times 
(t ~ 2t e ) where one expects the theory to break down 




Time 



FIG. 2: A comparison of simulated (unconnected dots) and 
predicted (solid lines) values of G2i(i) as a function of t. The 
results in the top row correspond to wave vectors pairs (1,2) 
and (2, 1) from left to right, the middle row correspond to 
(1, 3) and (3, 1), and the bottom row to (2, 3) and (3, 2). 
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due to the instantaneous approximation of the coupling 
vertices. 



V. SUMMARY 

In this paper, a mode-coupling theory was presented in 
which multiple point and multiple-time correlation func- 
tions for collective densities and tagged particle densities 
are expressed in terms of ordinary two-point, single-time 
interval correlation functions and a set of vertices. The 
theory developed here does not assume that fluctuating 
forces (noise) are Gaussian distributed and, in princi- 
ple, does not require an ansatz to obtain self-consistent 
equations. Furthermore, unlike kinetic theories, it is not 
restricted to low densities and should be applicable to 
dense fluids where cooperative motions of particles and 
collective modes are important. 

The formalism is based on projection operator tech- 
niques, which, for ordinary two-point, single-time corre- 
lation functions, lead to a generalized Langevin equation 
in which the memory function decays on a microscopic 
time scale (Eq.|SJ). The simple extension of the projection 
operator formalism to multiple-time correlation functions 
of tagged particle densities is complicated by the fact 
that the fluctuating forces appearing in the generalized 
Langevin equation do not obey Gaussian statistics. Fur- 
thermore, multiple-time correlations of the fluctuating 
force can in fact have a slow decay when the time ar- 
guments of these forces become comparable. In order 
to treat multiple-time correlation functions of fluctuat- 
ing forces properly, the correlation functions were ma- 
nipulated so that the time arguments of all fluctuating 
forces appearing in the correlations were guaranteed to 
be well-separated, ensuring that all memory functions 
which arise in the mode-coupling theory decay to zero on 
a molecular time scale. This construction allows equa- 
tions which are local in time to be obtained which relate 
the multiple-time correlation function to two-time but 
multiple-point correlations coupled by essentially time- 
independent vertices. These expressions, in turn, can be 
written as convolutions of two-point, single time inter- 
val "propagators" coupled by time- independent vertices. 
These propagators can either be taken directly from ex- 
periment, simulation, or can be solved self-consistently 
within the mode-coupling formalism. The vertices, which 
are composed of a static part (Euler term) and a gener- 
alized transport coefficient, can similarly be calculated 
from kinetic theory or taken from molecular dynamics 
and Monte-Carlo simulations. 

The equations for higher-order correlation functions 
contain an infinite sum of terms which can be made 
tractable for systems with a finite correlation length by 
applying a cumulant expansion technique pioneered by 
Oppenheim and co-workers ^| l2ll termed the N- 
ordering method. The method was applied to obtain the 
leading order and first mode-coupling corrections of ex- 
pressions for tagged-particle density multiple-time and 



mixed tagged/collective particle density multiple-point 
correlation functions designed to probe detailed aspects 
of the relaxation profile of glassy systems. 

The mode-coupling theory outlined here is the first 
step towards a fully microscopic theory applicable to 
molecular length scales that enables one to analyze the 
mapping of the dynamics in deterministic glassy systems 
to stochastic spin model of glasses, as well as to define 
sub-ensembles via filters based on dynamical and spatial 
properties to probe dynamic heterogeneities, clustering 
properties and many other aspects of the dynamics of 
frustrated systems. 
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APPENDIX: EVALUATION OF THE 
MODE-COUPLING VERTICES 

As mentioned in the main text, the vertices required 
for the calculation of the mode-coupling correction to 
the factorization of the multiple-time correlation function 
F 3 (h,t 2 ) are M^ Pl ' N \ M*** 1 and M^ H ' N \ where 
the general form of the vertex M|f is 

MH = M ii s + M 2i ss 
M^" '= W2QT) ■ Ktr 1 (AT) 

M£ ss = - dr (<p s 2 (T)<p?) ■ K{[-\ (A.2) 
Jo 

where is referred to as the "Euler" contribution to 
the vertex, M 2l is the "dissipative" part of the vertex, 
and ip(t) is the fluctuating force defined in Eq. EH 

Examining first the Euler contributions M 21 ss , we note 
that since {CQ 2 Q{*) = - (Q 2 CQf), we have 

M 2 t s (A> 1 (k-q)B(q);A> 1 (k)) = 

-(g^(A> 1 (k-q),B(q))Pt(k)\-k. (A.3) 

^From the time-reversal symmetry properties of the cor- 
relation functions, we see that M 21 ss — for B — N and 
B = H, and for B = Pi we find 

M 2 t s (A\(k- q)p,(q);^i(k)) = (A.4) 
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Evaluation of the dissipative vertices M^ ss requires 
the calculation of the linear ipi and bi-lincar ip 2 fluctu- 
ating forces. For the linear fluctuating force, we have (in 
the thermodynamic limit) 



iPi(k) • k/m, 



(A.5) 



and the bi-linear fluctuating force is 



^(k-q,q;t) = e F ^£^(k-q,q), 

where T' ± is a projection operator that projects onto a 
subspace orthogonal to that spanned by the multi-linear 
basis set Q a . Noting that (in the thermodynamic limit) 



Q s 2 (Ni(k - q)JV(q)) = 7Vi(k - q)iV(q) 

Q5(M(k-q)fl(q)) =JV 1 (k-q)fl(q) 
Q^N^k - q)ff(q)) = JVi(k - q)ff(q) 



g(q) 

(TV) 



JVi(k) 



(A.6) 



and "P^iV^k - q)P(q) = 0, we see that at i = 0, the 
bilinear fluctuating forces are 



^(TVi(k-q)iV(q)) 

^(JVi(k-q)fJ(q)) 
^(^(k-q^q)) 



i(k-q) 



m 



Pi(k-q)JV(q) 



5(q) ik 

(iV) m 



Pi(k) 



i(k - q) 



m 



•P ± P 1 (k-q)P(q)-q 



+ 7> ± iV 1 (k-q)/;P(q).q 

fc^.^ Pl (k-q)i/(q) 

m 

+ 7 T > ± A> 1 (k-q)/lff(q). 



Inserting these results into the expressions for the dissi- 
pative vertices and expanding the resulting expressions 
in powers of the base wave vector fco, one finds 



M. 



ssN 1 N;N 1 
21 



(k - q) • k 



Jo 



dT{P^(T)PfN 



M. 



ssN 1 P i ;N 1 
21 



M. 



ssN 1 H;N 1 _ 



, fc 2 s(q) 

m 2 (N) J 

O(k 3 ) 

(k - q) • k 



dr(PUr)Pi)+0(kl) 



■ k f 

- dr<Pf(r)PfJ0 
Jo 

— rdr(jH(r)Pi)+O(fc 3 ), 
n Jo 



rn 
k_q 

m 



where we have defined iq- jjj(q) = -2/3/V6 "P^zCi^q). 
Noting that the Green-Kubo expression for the self- 
diffusion coefficient m 2 D = J °° dt (Pf (t)Pf) and that 

(AZV) = d(A)/d(0n) and = -d(A)/d0, where /x is 
the chemical potential of the system, allows us to write 



, . f dD\ 5(0) 



L>fc 2 + O(fcg) 



ssNtH-.Nt _ (k — q) • k 



21 



V6 



dj3 J 



-— / rfr0iv(r)Pi)+O(fc 3 ). 



The dissipative vertex M^ NlH ' Nl contains a new trans- 
port coefficient J °° dr (jff( T )Pi) which could, in princi- 
ple, be evaluated by doing a short-time expansion or us- 
ing uncorrelated collisions in a kinetic theory. Nonethe- 
less, we expect this contribution to be very small since 
(^h(t)P) is strictly zero due to the projection operator 
. In the numerical work in the main text, we have set 
this term to zero. 

Using the Enskog form for the self-diffusion coefficient 
and the Carnahan-Starling equation of state, we find that 



3D \ 
8D\ 



D fj(2fj - 5) 5(0) 
g(a) 2(1 - f})* (N) 

D D 77(277 — 5) 35(0) 
'2p + g(a)2{l-n) i 2(3(N)- 



(A.7) 



Thus, the M 2 i vertices used in the main text are: 



J . /r ssN 1 N;N 1 _ DS(0) 



k-(k-q)r7(2r7-5) 2 



g(a) 2(1 - f)Y 



M%l(NiPf,Ni 



ik ■ q 



M. 



ssN 1 h-,n 1 _ £>(k - q) ■ k 



21 
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